setwd("C:/Users/tim/Documents/Distance Covariance Project/distancecovariance/R Code/")
#setwd("C:/Documents and Settings/ORIE user1/My Documents/Distance Covariance Project/R Code/")
#setwd("C:/Documents and Settings/1083117/My Documents/Distance Covariance Project/R Code")
#optional to get the analysis running...
source("Functions/source.all.R")
source.all() #this will produce errors if there are any uninstalled packages
install.all.packages() #install any missing packages
source.all() #source again...

###################################
#used only for model IV...
#generate the best models
#best.modelIIIs=NULL
#for (a in c(2,3)){
#	for(b in c(4,5,6)){	
#		for (c in c(7,8,9,10)){
#			best.modelIIIs=c(best.modelIIIs,paste("1",a,b,c,sep=","))
#		}	
#	}
#}
N=1000
n=100
lambda0=.7
distributions=c("normal",letters[1:18])
for(distribution in distributions){
	title=paste("Model Selections for " ,distribution," data, ","n=",n,", N=",N,sep="")
	best.model.chosen.b2=0
	best.model.chosen.dcov=0
	best.model.chosen.rank.dcov=0
	good.model.chosen.b2=0
	good.model.chosen.dcov=0
	good.model.chosen.rank.dcov=0
	moderate.model.chosen.b2=0
	moderate.model.chosen.dcov=0
	moderate.model.chosen.rank.dcov=0
	bad.model.chosen.b2=0
	bad.model.chosen.dcov=0
	bad.model.chosen.rank.dcov=0

	for(h in 1:N){
		z=NULL

		if(distribution=="normal"){
			for (i in 1:10){
				z=cbind(z,rnorm(n))
			}	
		}else{ #jordan distribution
			for(i in 1:10){
				z=cbind(z,rjordan(distribution,n))
			}
		}
			modelI   = cbind(z[,1],z[,2],z[,3],z[,1]+.5*z[,4],z[,2]+.7*z[,5],z[,3]+z[,6])
			modelII  = cbind(z[,1],z[,2],z[,3],z[,1]+.5*z[,4],z[,2]+.7*z[,5],z[,3]+0.5*z[,6])
			modelIII = cbind(z[,1],z[,2],z[,3],z[,1]+0.8*z[,2]+0.6*z[,4],z[,2]+0.7*z[,5],z[,3]+.5*z[,6])
			modelIV  = cbind(z[,1],z[,2],z[,2]+z[,3],z[,4],z[,4]+.75*z[,5],2*z[,4]+.75*z[,5]+1.5*z[,6],z[,7],z[,7]+.5*z[,8],2*z[,7]+.5*z[,8]+z[,9],3*z[,7]+z[,8]+z[,9]+z[,10])	
			
			

		#screeplot(model.princomp)
		#loadings(model.princomp)

		#modelI.princomp$sdev
		model = modelIII
		#employ method B2
		model.princomp = prcomp(model)
		
		lambdas=model.princomp$sdev
		loadings = as.matrix(model.princomp$rotation)

		k.minus.p.lambdas=as.numeric(which(lambdas<lambda0))
		k.minus.p.lambdas=k.minus.p.lambdas[rank(lambdas[k.minus.p.lambdas])]
		k.minus.p.vars = NULL
		for(j in k.minus.p.lambdas){
			#variable.to.reject = which.max(abs(loadings[,j]))
			variable.to.reject = which.max(loadings[,j])

			if (is.na(match(variable.to.reject,k.minus.p.vars))){
				#haven't already included this variable, so include it

			}else{
				#find the next biggest eigenvalue by taking subsets of the eigenvector column
				k=0
				loadings.column.subset=loadings[,j]
				while((k < dim(loadings)[1]) && (!is.na(match(variable.to.reject,k.minus.p.vars)))){
					#zero this variable loading
					loadings.column.subset[variable.to.reject]=0
					#variable.to.reject = which.max(abs(loadings.column.subset))
					variable.to.reject = which.max(loadings.column.subset)
					k=k+1
				}
			}
			k.minus.p.vars=c(k.minus.p.vars,variable.to.reject)
		}	

		p=as.numeric(which(lambdas>=lambda0))
		#take the variables not discarded
		p.vars=seq(1,dim(loadings)[1])[-k.minus.p.vars]


		Rms=NULL
		regressors = ""
		for (j in p.vars){ 
			variate = paste("model[",j,",]",sep="")
			regressors = ifelse(j==p.vars[length(p.vars)],paste(regressors, variate, sep = ""),paste(regressors, variate," + ", sep = ""))
		}
		for (j in k.minus.p.vars){
			formula.string = paste("model[", j,",]", " ~ ",sep="")		
			formula.string = as.character(paste(formula.string,regressors ,sep=""))
			regression.formula = formula(formula.string)
			Rms=c(Rms,sqrt(summary(lm(regression.formula))$r.sq))		
		}
		p.vars=p.vars[rank(p.vars)]
		#apply both dcov methods
		dcov.vars=get.dcov.vars(model,p.vars,k.minus.p.vars,rank=FALSE)$dcov.vars
		rank.dcov.vars=get.dcov.vars(model,p.vars,k.minus.p.vars,rank=TRUE)$dcov.vars	

		
		
		p.vars.string=""
		for(i in 1:length(p.vars)){
			p.vars.string=ifelse(i==length(p.vars),paste(p.vars.string,p.vars[i],sep=""),paste(p.vars.string,p.vars[i],",",sep=""))
		}

		dcov.vars.string=""
			for(i in 1:length(dcov.vars)){
				dcov.vars.string=ifelse(i==length(dcov.vars),paste(dcov.vars.string,dcov.vars[i],sep=""),paste(dcov.vars.string,dcov.vars[i],",",sep=""))
			}
		rank.dcov.vars.string=""
			for(i in 1:length(rank.dcov.vars)){
				rank.dcov.vars.string=ifelse(i==length(rank.dcov.vars),paste(rank.dcov.vars.string,rank.dcov.vars[i],sep=""),paste(rank.dcov.vars.string,rank.dcov.vars[i],",",sep=""))
			}

		best.modelIII.1="1,2,3"
		best.modelIII.2="1,2,6"
		best.modelIIIs=c(best.modelIII.1,best.modelIII.2)
		
		#actually used for model IV
		
		good.modelIIIs=c("1,5,6","1,3,5","2,4,6","2,3,4","3,4,5","4,5,6")
		#good.modelIIIs=c(" ")
		moderate.modelIIIs=c("1,3,4","1,4,6")
		#moderate.modelIIIs=c(" ")
		
		if(!is.na(match(p.vars.string,best.modelIIIs))){
			best.b2.increment = 1
		}else{
			best.b2.increment = 0
		}
		if(!is.na(match(dcov.vars.string,best.modelIIIs))){
			best.dcov.increment = 1
		}else{
			best.dcov.increment = 0
		}
		if(!is.na(match(rank.dcov.vars.string,best.modelIIIs))){
			best.rank.dcov.increment = 1
		}else{
			best.rank.dcov.increment = 0
		}
		if(!is.na(match(p.vars.string,good.modelIIIs))){
			good.b2.increment = 1
		}else{
			good.b2.increment = 0
		}
		if(!is.na(match(dcov.vars.string,good.modelIIIs))){
			good.dcov.increment = 1
		}else{
			good.dcov.increment = 0
		}
		if(!is.na(match(rank.dcov.vars.string,good.modelIIIs))){
			good.rank.dcov.increment = 1
		}else{
			good.rank.dcov.increment = 0
		}
		if(!is.na(match(p.vars.string,moderate.modelIIIs))){
			moderate.b2.increment = 1
		}else{
			moderate.b2.increment = 0
		}
		if(!is.na(match(dcov.vars.string,moderate.modelIIIs))){
			moderate.dcov.increment = 1
		}else{
			moderate.dcov.increment = 0
		}
		if(!is.na(match(rank.dcov.vars.string,moderate.modelIIIs))){
			moderate.rank.dcov.increment = 1
		}else{
			moderate.rank.dcov.increment = 0
		}

		bad.b2.increment = ifelse((best.b2.increment+good.b2.increment+moderate.b2.increment)==0,1,0);
		bad.dcov.increment = ifelse((best.dcov.increment+good.dcov.increment+moderate.dcov.increment)==0,1,0);
		bad.rank.dcov.increment = ifelse((best.rank.dcov.increment+good.rank.dcov.increment+moderate.rank.dcov.increment)==0,1,0);

		#aggregate results
		best.model.chosen.b2 = 	best.model.chosen.b2 + best.b2.increment
		best.model.chosen.dcov = best.model.chosen.dcov + best.dcov.increment
		best.model.chosen.rank.dcov = best.model.chosen.rank.dcov + best.rank.dcov.increment

		good.model.chosen.b2 = good.model.chosen.b2 + good.b2.increment
		good.model.chosen.dcov = good.model.chosen.dcov + good.dcov.increment
		good.model.chosen.rank.dcov = good.model.chosen.rank.dcov + good.rank.dcov.increment

		moderate.model.chosen.b2 = moderate.model.chosen.b2 + moderate.b2.increment
		moderate.model.chosen.dcov = moderate.model.chosen.dcov + moderate.dcov.increment
		moderate.model.chosen.rank.dcov = moderate.model.chosen.rank.dcov + moderate.rank.dcov.increment

		bad.model.chosen.b2 = bad.model.chosen.b2 + bad.b2.increment
		bad.model.chosen.dcov = bad.model.chosen.dcov + bad.dcov.increment
		bad.model.chosen.rank.dcov = bad.model.chosen.rank.dcov + bad.rank.dcov.increment
		trial.results.dataframe=as.data.frame(rbind(cbind("Best",h,best.b2.increment,best.dcov.increment,best.rank.dcov.increment),
							cbind("Good",h,good.b2.increment,good.dcov.increment,good.rank.dcov.increment),
							cbind("Moderate",h,moderate.b2.increment,moderate.dcov.increment,moderate.rank.dcov.increment),
							cbind("Bad",h,bad.b2.increment,bad.dcov.increment,bad.rank.dcov.increment)))
		names(trial.results.dataframe)=c("Subset","Trial","B2","B2+DCOV","B2+Rank DCOV")
							

		if(h==1){
			total.trial.results.dataframe=trial.results.dataframe
		}else{
			total.trial.results.dataframe=rbind(total.trial.results.dataframe,trial.results.dataframe)
		}
			
		
		#progress
		print(h/N)

	}
	
	barplot.dataframe = as.data.frame(rbind(cbind(best.model.chosen.b2,good.model.chosen.b2,moderate.model.chosen.b2,bad.model.chosen.b2),cbind(best.model.chosen.dcov,good.model.chosen.dcov,moderate.model.chosen.dcov,bad.model.chosen.dcov),cbind(best.model.chosen.rank.dcov,good.model.chosen.rank.dcov,moderate.model.chosen.rank.dcov,bad.model.chosen.rank.dcov)))
	names(barplot.dataframe)=c("Best","Good","Moderate","Bad")
	dev.new()
	barplot(as.matrix(barplot.dataframe), main=title, ylab= "Total",   beside=TRUE, col=c("dark gray","gray","light gray"), legend.text=c("B2","B2+DCOV","B2+Rank DCOV"))
	write.csv(x=barplot.dataframe,file=paste(distribution,"_aggregate.csv",sep=""),append=FALSE,row.names = FALSE,col.names=TRUE)
	write.table(total.trial.results.dataframe,file=paste(distribution,"_trials.csv",sep=""),sep=",",append=FALSE,row.names = FALSE)
		
}